﻿#include "datasourcegdaldem.h"
#include "mainwindow.h"

#include <gdal_priv.h>
#include <QDebug>
#include <osgDB/WriteFile>

#include <qmath.h>

#include "GeoRaster.h"
using namespace rxtiler;
rxtiler::DataSourceGDALDEM::DataSourceGDALDEM()
{


}

rxtiler::DataSourceGDALDEM::~DataSourceGDALDEM()
{
    this->close();
}

bool DataSourceGDALDEM::open()
{
    _xResolution=0.0;
    _yResolution=0.0;
//    GDALAllRegister();//注册驱动
//    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");//支持中文路径
    poDataset = ( GDALDataset* )GDALOpen( _path.toStdString().c_str(), GA_ReadOnly );

    double adfGeoTransform[6];
    CPLErr res = poDataset->GetGeoTransform(adfGeoTransform);
    if (res == CE_None)
    {
        _xResolution=adfGeoTransform[1];
        _yResolution=-adfGeoTransform[5];
    }

    int src_bandsize = poDataset->GetRasterCount();
    if(src_bandsize < 1){
                GDALClose(poDataset);
                poDataset = 0;
                return false;
    }
    return true;
}

bool DataSourceGDALDEM::close()
{
    if(poDataset){
        GDALClose(  poDataset );
    }
    poDataset = 0;
    return true;
}



rxtiler::Extent rxtiler::DataSourceGDALDEM::getExtent()
{
    double adfGeoTransform[6];
    if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None){
        int x=poDataset->GetRasterXSize();
        int y=poDataset->GetRasterYSize();
        double left=adfGeoTransform[0];
        double top=adfGeoTransform[3];
        double right=adfGeoTransform[0]+x*adfGeoTransform[1]+y*adfGeoTransform[2];
        double bottom=adfGeoTransform[3]+x*adfGeoTransform[4]+y*adfGeoTransform[5];
        _tif.left = left;
        _tif.right = right;
        _tif.bottom = bottom;
        _tif.top = top;

        return rxtiler::Extent(left,right,bottom,top);
    }else
    {
        return rxtiler::Extent(0,0,0,0);
    }
}

bool rxtiler::DataSourceGDALDEM::computerLevels(rxtiler::TileOption *option, int *min, int* max)
{

    return true;
}

osg::HeightField *DataSourceGDALDEM::getDEM(TileKey &key, Extent &extent)
{
    int dem_size = TileOption::outTileOption()->tilesize;
    if(!_tif.cover(extent)){
        return 0;
    }

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();

    Extent ext2 = extent.intersect(_tif);

    int x = (iWidth-1) * (ext2.left - _tif.left) / _tif.width() - 1;
    int y = (iHeight-1) * ( _tif.top - ext2.top) / _tif.height() - 1;
    int XSize = ::ceil((iWidth-1) * ext2.width() /  _tif.width()) + 2;
    int YSize = ::ceil((iHeight-1) * ext2.height()  /  _tif.height()) + 2;


    if(x < 0) x = 0;
    if(y < 0) y = 0;
    if(x + XSize >= iWidth) XSize = XSize - 1;
    if(y + YSize >= iHeight) YSize = YSize -1;

    if(XSize < 3 || YSize< 3) return 0;

    ext2.left = _tif.left + _tif.width() * x / (iWidth-1);
    ext2.right = _tif.left + _tif.width() * (x+XSize) / (iWidth-1);
    ext2.top = _tif.top - _tif.height() * y / (iHeight-1);
    ext2.bottom = _tif.top- _tif.height() * (y+YSize) / (iHeight-1);


    int bxSize = XSize;
    int bySize = YSize;
    if(bxSize > 256) bxSize = 256;
    if(bySize > 256) bySize = 256;

    float* buffer = (float*)malloc(bxSize * bySize * 4);
    poDataset->RasterIO(GF_Read,x,y,XSize,YSize,buffer,bxSize,bySize,GDT_Float32,1,0,0,0,0);

    osg::ref_ptr<osg::HeightField> dem1 = new osg::HeightField();
    dem1->allocate(bxSize,bySize);
    for(int x =0; x< bxSize;x++) {
        for(int y = 0; y< bySize; y++){
            float h=buffer[bxSize*(bySize -y -1)+x];
            dem1->setHeight(x,y,h);
        }
    }
    free(buffer);
    if(ext2.equal(extent) && bxSize == dem_size && bySize == dem_size){
        return dem1.release();
    }
    osg::HeightField* dem2 = new osg::HeightField();
    dem2->allocate(dem_size,dem_size);
    mapRaster(dem1,ext2,dem2,extent);

    return dem2;
}

osg::Image *rxtiler::DataSourceGDALDEM::getImage(TileKey &key, Extent &extent)
{
    return 0;
}

/*
osg::HeightField *rxtiler::DataSourceGDAL::getDEM(rxtiler::TileKey &key, rxtiler::Extent &extent)
{
    int dem_size = TileOption::outTileOption()->tilesize;
    if(!_tif.cover(extent)){
        return 0;
    }

    int iWidth = poDataset->GetRasterXSize();
    int iHeight = poDataset->GetRasterYSize();

    Extent ext2 = extent.intersect(_tif);

    int x = (iWidth-1) * (ext2.left - _tif.left) / _tif.width() - 1;
    int y = (iHeight-1) * ( _tif.top - ext2.top) / _tif.height() - 1;
    int XSize = ::ceil((iWidth-1) * ext2.width() /  _tif.width()) + 2;
    int YSize = ::ceil((iHeight-1) * ext2.height()  /  _tif.height()) + 2;


    if(x < 0) x = 0;
    if(y < 0) y = 0;
    if(x + XSize >= iWidth) XSize = XSize - 1;
    if(y + YSize >= iHeight) YSize = YSize -1;

    if(XSize < 3 || YSize< 3) return 0;

    ext2.left = _tif.left + _tif.width() * x / (iWidth-1);
    ext2.right = _tif.left + _tif.width() * (x+XSize) / (iWidth-1);
    ext2.top = _tif.top - _tif.height() * y / (iHeight-1);
    ext2.bottom = _tif.top- _tif.height() * (y+YSize) / (iHeight-1);


    int bxSize = XSize;
    int bySize = YSize;
    if(bxSize > 1024) bxSize = 1024;
    if(bySize > 1024) bySize = 1024;

    float* buffer = (float*)malloc(bxSize * bySize * 4);
    poDataset->RasterIO(GF_Read,x,y,XSize,YSize,buffer,bxSize,bySize,GDT_Float32,1,0,0,0,0);

    osg::ref_ptr<osg::HeightField> dem1 = new osg::HeightField();
    dem1->allocate(bxSize,bySize);
    for(int x =0; x< bxSize;x++) {
        for(int y = 0; y< bySize; y++){
            float h=buffer[bxSize*(bySize -y -1)+x];
            dem1->setHeight(x,y,h);
        }
    }
    free(buffer);
    if(ext2.equal(extent) && bxSize == dem_size && bySize == dem_size){
        return dem1.release();
    }
    osg::HeightField* dem2 = new osg::HeightField();
    dem2->allocate(dem_size,dem_size);
    mapRaster(dem1,ext2,dem2,extent);

    return dem2;

}
*/
